Precession Of The Equinoxes

1. 外力矩

要推导岁差, 第一步需要得到地球所受的外力矩, 而要推导地球所受引力矩的表达式, 需分析太阳或月球对地球赤道隆起的引力作用.

假设地球是均质扁球体 (赤道转动惯量 A, 自转轴转动惯量 C > A), 外引力源 (质量 M, 距离 d) 的方向 $\hat{r}$ 与地球自转轴 $\hat{z}$ 的夹角为黄赤交角 $ε \approx {23.44}^{\circ}$ .

地球转动惯量张量在自转系中写成对角矩阵:

I = (\begin{matrix} A & 0 & 0 \\ 0 & A & 0 \\ 0 & 0 & C \end{matrix}), Δ I = C - A > 0

1.1 引力势展开

引力源在地球内部位置 ${\vec{r}}^{'}$ 处产生的引力势是:

Φ ({\vec{r}}^{'}) = - \frac{G M}{| \vec{d} - {\vec{r}}^{'} |}, \vec{d} = d \hat{r}

在 | ${\vec{r}}^{'} | ≪ d$ 的近似下, 将上式展开至二阶项:

Φ ({\vec{r}}^{'}) \approx - \frac{G M}{d} [1 + \frac{{\vec{r}}^{'} \cdot \hat{r}}{d} + \frac{3 (\hat{r} \cdot {\vec{r}}^{'})^{2} - r^{' 2}}{2 d^{2}}]

其中常数项和偶极项对净力矩无贡献, 主导项为二阶项:

Φ_{2} ({\vec{r}}^{'}) = - \frac{G M}{2 d^{3}} [3 (\hat{r} \cdot {\vec{r}}^{'})^{2} - r^{' 2}]

1.2 引力展开和微元受力

已知引力场 $\vec{g} = - \nabla^{'} Φ$ , 质量元 $d m$ 受力 $d \vec{F} = - \nabla^{'} Φ_{2} d m$ , $\nabla^{'}$ 表示关于 ${\vec{r}}^{'}$ 求梯度, 可以得到:

\nabla [3 (\hat{r} \cdot {\vec{r}}^{'})^{2} - r^{' 2}] = 6 (\hat{r} \cdot {\vec{r}}^{'}) \hat{r} - 2 {\vec{r}}^{'}

因此:

d \vec{F} = \frac{G M}{d^{3}} [3 (\hat{r} \cdot {\vec{r}}^{'}) \hat{r} - {\vec{r}}^{'}] d m

1.3 总力矩积分

积分得到地球所受总力矩 $\vec{τ} = \int {\vec{r}}^{'} \times d \vec{F}$ :

\vec{τ} = \frac{G M}{d^{3}} \int {\vec{r}}^{'} \times [3 (\hat{r} \cdot {\vec{r}}^{'}) \hat{r} - {\vec{r}}^{'}] d m = \frac{3 G M}{d^{3}} \int ({\vec{r}}^{'} \times \hat{r}) (\hat{r} \cdot {\vec{r}}^{'}) d m

其中:

\int (\hat{r} \cdot {\vec{r}}^{'}) {\vec{r}}^{'} d m = I \cdot \hat{r} = (\begin{matrix} A r_{x} \\ A r_{y} \\ C r_{z} \end{matrix})

r_{z} = \hat{r} \cdot \hat{z} = \sin ε

于是原积分可以改写为:

\begin{aligned} \int (\hat{r} \cdot {\vec{r}}^{'}) ({\vec{r}}^{'} \times \hat{r}) d m & = (\int (\hat{r} \cdot {\vec{r}}^{'}) {\vec{r}}^{'} d m) \times \hat{r} \\ = (\begin{array}{c} A r_{x} \\ A r_{y} \\ C r_{z} \end{array}) \times (\begin{array}{c} r_{x} \\ r_{y} \\ r_{z} \end{array}) = (C - A) r_{z} (\begin{array}{c} - r_{y} \\ r_{x} \\ 0 \end{array}) \end{aligned}

其中 $(\begin{matrix} - r_{y} \\ r_{x} \\ 0 \end{matrix}) = \hat{z} \times \hat{r}$ , 所以:

\int (\hat{r} \cdot {\vec{r}}^{'}) ({\vec{r}}^{'} \times \hat{r}) d m = (C - A) \sin ε (\hat{z} \times \hat{r})

最终得到:

\vec{τ} = \frac{3 G M}{d^{3}} (C - A) \sin ε (\hat{z} \times \hat{r})

τ = \frac{3 G M}{2 d^{3}} (C - A) \sin 2 ε

但是, 天体在黄道上是运动的, 这使得我们直接使用 $ε$ 这样一个定值是不准确的. 实际上, 将 $ε$ 视为天体的瞬时赤经, 很容易将上式改写成:

\vec{τ} = - \frac{3 G M}{d^{3}} (C - A) \sin δ \cos δ (\begin{matrix} \sin α \\ - \cos α \\ 0 \end{matrix})

其中的 $δ, α$ 是天体瞬时的赤纬和赤经, 取值范围为 $δ \in [- ε, ε], α \in [- 180^{\circ}, 180^{\circ})$ . 上式中, $\sin δ \cos α$ 明显在长时间平均中为 0, 而 $\sin δ \sin α$ 始终保持不变的符号(正或负), 根据球面直角三角形的规律, 可以发现:

\begin{array}{r} \sin δ = \sin ε \sin u \\ \cos δ = \cos u / c o s α \\ \tan α = \cos ε \tan u \end{array}} ⟹ {\begin{array}{r} \sin δ \cos δ \sin α = \sin ε \cos ε \sin^{2} u \\ \sin δ \cos δ \cos α = \sin ε \sin u \cos u \end{array}

u 是天体当前的春分点角距(沿着黄道到春分点的角度, 或者叫黄经, 如图), 它正比于时间 t, 可以发现 x 方向的力矩在长时间平均后为 $\frac{1}{2}$ , 即:

\frac{1}{2 π} \int_{0}^{2 π} \sin^{2} u d u = \frac{1}{2}

所以, 最终得到的力矩需要加上 $\frac{1}{2}$ 倍的修正, 最终得到长时间平均后的力矩是:

⟨ \vec{τ} ⟩ = \frac{3 G M}{2 d^{3}} (C - A) \cos ε \sin ε

2. 进动角速度

2.1 角动量进动和外力力矩的关系

如下图, 不需要考虑 $z$ 方向力矩, 它只对物体的自转起作用. 所以只考虑垂直于 $z$ 方向的力矩 $M$ .

考虑 x 方向的力矩 M, 根据力矩等于角动量变化率得到:

Δ J = J s i n θ Δ φ = M Δ t \to Ω = \frac{Δ φ}{Δ t} = \frac{M}{J \sin θ}

可以发现 x 方向的力矩带来的效应是角动量以角速度 $Ω$ 进行进动. 类似的, y 方向的力矩将使得角动量发生章动.

2.2 岁差

结合前面的结果, 有:

\begin{aligned} Ω & = \frac{τ}{J \sin ε} \\ τ & = \frac{3 G M}{4 d^{3}} (C - A) \sin 2 ε \\ J & \approx C ω \end{aligned}

$ω$ 为地球的自转角速度. 也即:

Ω \approx \frac{3 G M}{d^{3}} \frac{C - A}{2 C ω} \cos ε = 3 G \frac{C - A}{2 C ω} (\frac{M_{⊙}}{d_{⊙}^{3}} \cos ε_{⊙} + \frac{M_{☾}}{d_{☾}^{3}} \cos ε_{☾})

通过计算(见[[#用 ipython 计算结果的代码]])得到它的结果是 50.45''/a, 对应的周期是两万五千七百年, 与天文观测的结果接近.

参考

https://en.m.wikipedia.org/wiki/Axial_precession

用 ipython 计算结果的代码

import math

# 物理常数
G = 6.674e-11  # 万有引力常数 (m^3 kg^-1 s^-2)
omega = 7.272e-5  # 地球自转角速度 (rad/s)
epsilon_deg_sun = 23.4392  # 黄赤交角 (degrees)
epsilon_rad_sun = math.radians(epsilon_deg_sun)  # 黄赤交角 (radians)
cos_epsilon_sun = math.cos(epsilon_rad_sun)
epsilon_deg_moon = epsilon_deg_sun + 5.1453  # 白赤交角
epsilon_rad_moon = math.radians(epsilon_deg_moon)  # 黄赤交角 (radians)
cos_epsilon_moon = math.cos(epsilon_rad_moon)

# 地球参数
C = 8.034e37  # 地球绕极轴转动惯量 (kg·m^2)
f = 1 / 298.257  # 地球扁率
C_minus_A = f * C  # C - A

# 太阳和月球参数
M_sun = 1.989e30  # 太阳质量 (kg)
M_moon = 7.342e22  # 月球质量 (kg)
d_sun = 1.496e11  # 日地平均距离 (m)
d_moon = 3.844e8  # 地月平均距离 (m)

# 计算太阳和月球的引力项
sun_term = M_sun / (d_sun ** 3) * cos_epsilon_sun
moon_term = M_moon / (d_moon ** 3)* cos_epsilon_moon
sum_term = sun_term + moon_term

# 计算转动惯量相关系数
coefficient = (3 * G * C_minus_A) / (2 * C * omega)

# 计算岁差角速度 (rad/s)
omega_precession = coefficient * sum_term

# 单位转换
seconds_per_year = 3.154e7  # 一年的秒数
arcsec_per_rad = 206265  # 每弧度对应的角秒数
omega_precession_arcsec_per_year = omega_precession * seconds_per_year * arcsec_per_rad

print(f"太阳引力项: {sun_term:.2e} kg/m³")
print(f"月球引力项: {moon_term:.2e} kg/m³")
print(f"引力项总和: {sum_term:.2e} kg/m³")
print(f"转动惯量系数: {coefficient:.2e} m³/(kg·s)")
print(f"岁差角速度: {omega_precession:.2e} rad/s")
print(f"岁差速率: {omega_precession_arcsec_per_year:.2f} 角秒/年")
print(f"岁差: {360 * 60 * 60 / omega_precession_arcsec_per_year:.2f} 年")
   
--------------------------------
   
太阳引力项: 5.45e-04 kg/m³
月球引力项: 1.14e-03 kg/m³
引力项总和: 1.68e-03 kg/m³
转动惯量系数: 4.62e-09 m³/(kg·s)
岁差角速度: 7.75e-12 rad/s
岁差速率: 50.45 角秒/年
岁差: 25689.14 年